2. Magnetostatics

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import math

2.1. model of the coil:

cyl = Cylinder((0,0,0), Z, r=0.01, h=0.03).faces[0]
heli = Edge(Segment((0,0), (12*math.pi, 0.03)), cyl)
ps = heli.start
vs = heli.start_tangent
pe = heli.end
ve = heli.end_tangent

e1 = Segment((0,0,-0.03), (0,0,-0.01))
c1 = BezierCurve( [(0,0,-0.01), (0,0,0), ps-vs, ps])
e2 = Segment((0,0,0.04), (0,0,0.06))
c2 = BezierCurve( [pe, pe+ve, (0,0,0.03), (0,0,0.04)])
spiral = Wire([e1, c1, heli, c2, e2])
circ = Face(Wire([Circle((0,0,-0.03), Z, 0.001)]))
coil = Pipe(spiral, circ)

coil.faces.maxh=0.2
coil.faces.name="coilbnd"
coil.faces.Max(Z).name="in"
coil.faces.Min(Z).name="out"
coil.mat("coil")
crosssection = coil.faces.Max(Z).mass

Draw (coil);

box = Box((-0.04,-0.04,-0.03), (0.04,0.04,0.06))
box.faces.name = "outer"
air = box-coil
air.mat("air");

2.2. mesh-generation of coil and air-box:

geo = OCCGeometry(Glue([coil,air]))
with TaskManager():
    mesh = Mesh(geo.GenerateMesh(meshsize.coarse, maxh=0.01)).Curve(3)
    
Draw (mesh, clipping={"y":1, "z":0, "dist":0.012});

checking mesh data materials and boundaries:

mesh.ne, mesh.nv, mesh.GetMaterials(), mesh.GetBoundaries()

(153516,
 26451,
 ('coil', 'air'),
 ('out',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'coilbnd',
  'in',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer'))

2.3. Solve a potential problem to determine current density in wire:

on the domain \(\Omega_{\text{coil}}\):

\[\begin{eqnarray*} j & = & \sigma \nabla \Phi \\ \operatorname{div} j & = & 0 \end{eqnarray*}\]

port boundary conditions:

\[\begin{eqnarray*} \Phi & = & 0 \qquad \qquad \text{on } \Gamma_{\text{out}}, \\ j_n & = & \frac{1}{|S|} \quad \qquad \text{on } \Gamma_{\text{in}}, \end{eqnarray*}\]

and \(j_n=0\) else

fespot = H1(mesh, order=3, definedon="coil", dirichlet="out")
phi,psi = fespot.TnT()
sigma = 58.7e6
bfa = BilinearForm(sigma*grad(phi)*grad(psi)*dx).Assemble()
inv = bfa.mat.Inverse(freedofs=fespot.FreeDofs(), inverse="sparsecholesky")
lff = LinearForm(1/crosssection*psi*ds("in")).Assemble()
gfphi = GridFunction(fespot)
gfphi.vec.data = inv * lff.vec

Draw (gfphi, draw_vol=False, clipping={"y":1, "z":0, "dist":0.012});

2.4. Solve magnetostatic problem:

current source is current from potential equation:

\[ \int \mu^{-1} \operatorname{curl} u \cdot \operatorname{curl} v \, dx = \int j \cdot v \, dx \]

fes = HCurl(mesh, order=2, nograds=True)
print ("HCurl dofs:", fes.ndof)
u,v = fes.TnT()
mu = 4*math.pi*1e-7
a = BilinearForm(1/mu*curl(u)*curl(v)*dx+1e-6/mu*u*v*dx)
pre = Preconditioner(a, "bddc")
f = LinearForm(sigma*grad(gfphi)*v*dx("coil"))
with TaskManager():
    a.Assemble()
    f.Assemble()

HCurl dofs: 796070

from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates=True)
gfu = GridFunction(fes)
with TaskManager():
    gfu.vec.data = inv * f.vec

CG iteration 1, residual = 23.480063276896825     
CG iteration 2, residual = 0.08454523020404824     
CG iteration 3, residual = 0.017179588184440417     

CG iteration 4, residual = 0.009823452055571237     
CG iteration 5, residual = 0.006174780720300533     
CG iteration 6, residual = 0.004093419442077609     

CG iteration 7, residual = 0.0028077125091493875     
CG iteration 8, residual = 0.0017770290568212093     
CG iteration 9, residual = 0.0012198681828963344     

CG iteration 10, residual = 0.0009383538181344706     
CG iteration 11, residual = 0.000726160979221225     
CG iteration 12, residual = 0.0005721827343848372     

CG iteration 13, residual = 0.0004689027415680588     
CG iteration 14, residual = 0.0004101649169496657     
CG iteration 15, residual = 0.0003551390248698898     

CG iteration 16, residual = 0.00027787663143305437     
CG iteration 17, residual = 0.00021440113770878052     
CG iteration 18, residual = 0.00016101476248652627     

CG iteration 19, residual = 0.00011727050090439413     
CG iteration 20, residual = 8.225595406137369e-05     
CG iteration 21, residual = 5.6606813579813626e-05     

CG iteration 22, residual = 4.0639974746485745e-05     
CG iteration 23, residual = 2.683637731404489e-05     
CG iteration 24, residual = 1.895633593089709e-05     

CG iteration 25, residual = 1.3147181429820023e-05     
CG iteration 26, residual = 9.491050169650074e-06     
CG iteration 27, residual = 7.269688907783154e-06     

CG iteration 28, residual = 5.557837986141231e-06     
CG iteration 29, residual = 4.00338787718912e-06     
CG iteration 30, residual = 2.8614289191086017e-06     

CG iteration 31, residual = 2.2253960281128615e-06     
CG iteration 32, residual = 1.9497875704788573e-06     
CG iteration 33, residual = 1.4786105339822238e-06     

CG iteration 34, residual = 1.1955495929251882e-06     
CG iteration 35, residual = 9.310037425851402e-07     
CG iteration 36, residual = 7.755618490373793e-07     

CG iteration 37, residual = 6.064785767934925e-07     
CG iteration 38, residual = 4.6367240106339713e-07     
CG iteration 39, residual = 3.5480777833271324e-07     

CG iteration 40, residual = 2.577679525648949e-07     
CG iteration 41, residual = 1.84443210712618e-07     
CG iteration 42, residual = 1.3425502392336177e-07     

CG iteration 43, residual = 1.0034181783285406e-07     
CG iteration 44, residual = 7.424387151387077e-08     
CG iteration 45, residual = 5.4521421244939066e-08     

CG iteration 46, residual = 4.2187140208235934e-08     
CG iteration 47, residual = 3.382972513367436e-08     
CG iteration 48, residual = 2.7655500667671966e-08     

CG iteration 49, residual = 2.234252099708594e-08     
CG iteration 50, residual = 1.683141828801914e-08     
CG iteration 51, residual = 1.2614023654736312e-08     

CG iteration 52, residual = 1.0204803239724687e-08     
CG iteration 53, residual = 8.49138396960304e-09     
CG iteration 54, residual = 6.941385234162937e-09     

CG iteration 55, residual = 5.1849153334310215e-09     
CG iteration 56, residual = 3.911147900873985e-09     
CG iteration 57, residual = 3.024197598816761e-09     

CG iteration 58, residual = 2.3871447995872257e-09     
CG iteration 59, residual = 1.731155777449626e-09     
CG iteration 60, residual = 1.3070496018401746e-09     

CG iteration 61, residual = 9.263920184022558e-10     
CG iteration 62, residual = 6.861337941193581e-10     
CG iteration 63, residual = 5.344917595039201e-10     

CG iteration 64, residual = 4.219511731258286e-10     
CG iteration 65, residual = 3.25715465784038e-10     
CG iteration 66, residual = 2.5602275714343417e-10     

CG iteration 67, residual = 2.1564218256566896e-10     
CG iteration 68, residual = 1.8368248896806872e-10     
CG iteration 69, residual = 1.4713836107355411e-10     

CG iteration 70, residual = 1.1069363418470037e-10     
CG iteration 71, residual = 8.631587758762912e-11     
CG iteration 72, residual = 6.79248260594301e-11     

CG iteration 73, residual = 5.88387793706594e-11     
CG iteration 74, residual = 4.2020043715037714e-11     
CG iteration 75, residual = 3.069644750752258e-11     

CG iteration 76, residual = 2.3155758116110407e-11     

# Draw (curl(gfu), mesh, draw_surf=False, \
#      min=0, max=3e-4, clipping = { "y":1, "z" : 0, "function":False}, vectors = { "grid_size":100});

# from ngsolve import *
# from netgen.occ import unit_cube
# from ngsolve.webgui import Draw
# mesh = Mesh (unit_cube.GenerateMesh(maxh=0.2))
#
# cf = CF((-y,x,0.1))
N = 5
# p = [(i/N,0.1,k/N) for i in range(1,N) for k in range(1,N)]
p = [(-0.04*0.08*i/N,0.0,-0.04+0.08*k/N) for i in range(1,N) for k in range(1,N)]
fieldlines = curl(gfu)._BuildFieldLines(mesh, p, (N-1)**2)
Draw(curl(gfu), mesh, draw_vol=False, objects=[fieldlines], draw_surf=False, \
     min=0, max=3e-4, autoscale=False, clipping = { "y":1, "z" : 0, "function":False} )

Fieldline Calculation 0%
Fieldline Calculation 6%
Fieldline Calculation 12%
Fieldline Calculation 18%
Fieldline Calculation 25%
Fieldline Calculation 31%
Fieldline Calculation 37%
Fieldline Calculation 43%
Fieldline Calculation 50%
Fieldline Calculation 56%
Fieldline Calculation 62%
Fieldline Calculation 68%
Fieldline Calculation 75%
Fieldline Calculation 81%
Fieldline Calculation 87%
Fieldline Calculation 93%
Fieldline Calculation 100%

BaseWebGuiScene